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Abstract 



We perform cosmological perturbation theory in Hassan-Rosen bimetric gravity for general homo- 
geneous and isotropic backgrounds. In the de Sitter approximation, we obtain decoupled sets of 
massless and massive scalar gravitational fluctuations. Matter perturbations then evolve like in 
Einstein gravity. We perturb the future de Sitter regime by the ratio of matter to dark energy, 
producing quasi-de Sitter space. In this more general setting the massive and massless fluctuations 
mix. We argue that in the quasi-de Sitter regime, the growth of structure in bimetric gravity differs 
from that of Einstein gravity. 
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1 Introduction 

The cosmological constant problem [1] [2j is one of the most vexing problems in physics. The main 
problem is why the vacuum energy densities of quantum field theory seem to contribute to observable 
gravitational physics so much less than simple estimates would indicate. Presumably, the problem 
would be resolved in a theory of quantum gravity. String theory contains some quantum gravity, 
but it is notoriously difficult to address quantum problems in dynamical gravity with the present 
formulation of string theory. With the observation of the accelerated expansion of the universe in 
1998, usually attributed to a dominant dark energy component such as the cosmological constant (or 
quintessence, which has similar naturalness problems), the issue has been driven to a sharp point. 
As reviewed in [SJ, very few of the many proposed solutions to the problem stand a remote chance 
of success. 

Modified gravity is one of the well-known proposed "solutions" that fares particularly poorly 
in the evaluation of e.g. [2], because modified gravity theories are typically only deep-infrared 
modifications of gravity (where "deep-infrared" means very low energies, some tiny fraction of an 
electron volt), which ultimately seems insufficient to solve the problem of quantum field theory 
contributions from all known fields, including for example around the electron mass of 511 keV. 
There are suggestions how modified gravity could effectively limit how energy gravitates (a "filter" ) 
at a wider variety of energy scales, as in the proposed "degravitation" mechanism [3j, and the 
earlier discussions of screening mechanisms summarized in [4j. These are intriguing but incomplete 
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suggestions, in that it is not yet clear if any of these mechanisms are actually realized in any 
underlying theory in which the range of applicability of these mechanisms could be reliably evaluated. 

Massive gravity, a theory of gravity where the graviton has a mass (which is typically constrained 
by observations to be extremely small, perhaps 10 -33 eV, see e.g. [5] for a list of references) is 
at face value a relatively minor and again deep-infrared modification of gravity. The study of 
massive gravity was initiated in 1939 by Fierz and Pauli [7], but these theories suffered from ghost 
instabilities at the nonlinear level. In 2010 progress was made when a particular set of nonlinear 
ghost-free interactions was found by de Rham, Gabadadze and Tolley (dRGT) in a series of papers 
H2HgQ33|2iai2Il|221|231|23| , following seminal earlier work in [Ml [39]. The dRGT formulation 
of nonlinear massive gravity requires a fixed auxiliary two-tensor /„„ with no dynamics of its own. 
Apart from aesthetic concerns about this, for our purposes it is a deficiency of this theory that it 
seems to have no homogeneous and isotropic cosmological solutions [181 ED] • 

Last year, Hassan and Rosen |14|. [T5l [16] gave dynamics to the tensor f^ u in a bimetric theory, 
with the nonlinear interactions between g^ u and imported from the de Rham-Gabadadze- Tolley 
massive gravity theories. (In fact, the need for two metrics to realize massive gravity covariantly 
was already appreciated in 1976 [13], but there was no theory without nonlinear instabilities.) The 
Hassan- Rosen bimetric theory has cosmological solutions, as explored in [81 [91 [XT] HOj . In the 
aforementioned cosmological solutions of Hassan-Rosen bimetric gravity, both g^ v and f^ u have 
equations of motion, and the background solutions we consider are of general FLRW form for both 
backgrounds. 

Now, bimetric gravity is a more far-reaching modification of gravity than massive gravity, in that 
the new gravitationally coupled tensor field f^ u has dynamics of its own. (In the formulation we will 
be using here, it does not couple directly to matter, so "bimetric" is a little bit of a misnomer, "gravity 
coupled to matter and a symmetric two-tensor" would have been more accurate.) Of course, since 
the field content and interactions of bimetric theory are different from Einstein gravity, this theory 
may have different quantum properties at any given scale and not exclusively in the deep infrared. 
One way to try to understand the quantum properties of this theory would be to try to embed the 
theory in string theory, but currently it is not known how to do this. On the good side, the Hassan- 
Rosen bimetric theory (respectively de Rham-Gabadadze- Tolley massive gravity theory) has the 
kind of rigid structure that one would think could possibly descend from an underlying theory, like 
string theory. Symmetries constrain the interaction terms to the form V(f~ 1 g), and their relative 
coefficients are constrained, and it is now understood how to construct these theories in various 
dimensions and including higher-derivative corrections [29]. It would be somewhat surprising, and 
a shame, if this structure existed for no reason at all. 

To be clear, there is so far no clear indication that even embedding Hassan-Rosen bimetric 
theory in string theory would particularly help with the cosmological constant problem, but at 
least the problems could be addressed in a theory that is apparently nonlinearly consistent and also 
fundamentally different from Einstein gravity already at the level of the low-energy effective action. 

On the other hand, it may be easier to rule these kinds of theories out observationally (and 
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classically) than to properly understand their quantization, so here we pursue strategies to achieve 
classical observational tests. In this paper, we 



derive the linearized gravitational scalar fluctuation equations 
for FLRW backgrounds (section [4]) 

find convenient gauge invariant variables (section [4| 

solve the equations in the special case of a de Sitter background (section [6]) 
— this is not completely new, see [25j ES] 

develop the quasi-de Sitter (qdS) approximation in bimetric gravity (section [7]) 

find solutions of the qdS fluctuation equations, both analytical and numerical (section [7]) 

in general, construct some necessary framework for the analysis of growth of structure in 
Hassan-Rosen bimetric gravity. 



Detailed observational and phenomenological analyses are left for the future. 

We also mention that there has also been recent related work on multi-metric theory, the natural 
generalization of this framework to coupling multiple spin- two fields nonlinear ly [37} l35l [36] . 

Finally, there is also progress on related theories in three dimensions [33]. In fact, the "new 
massive gravity" theory in three dimensions [32J, which generated some excitement in the last few 
years, is a scaling limit of the Hassan- Rosen bimetric theory [34} 129]. 



2 Hassan-Rosen bimetric massive gravity 

The bimetric massive gravity theory found by Hassan and Rosen [15J is given by the action 

M 2 r . M 2 r — 

Sur = —± J d 4 x^gR(g) - -f J d 4 x^~fR(f) (l) 

4 

+ m 2 M 2 g I' d A xV=g P^n ( V*? 31 /) + / d A x^C m (g, <f) . (2) 

J n=0 ^ 

and represents a natural generalization of the de Rham-Gabadadze-Tolley massive gravity theory [22] 
to a theory with two dynamical metrics, as discussed in the introduction. Here (5 n are free parameters, 
which in general are the coefficients in the "deformed determinant" of [12] . The interaction terms 
e n (X) are elementary symmetric polynomials of the eigenvalues of the matrix X, which explicitly 
are given by 

e (X) = l, ei (X) = TrX, e 2 (X) = ± ((Tr X) 2 - Tr X 2 ) , 
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e 3 (X) = - ( (Tr X) a -3TVXTrI 2 + 2TrX 3 
6 



e 4 (X) = det (X) 



(3) 



We have chosen to only couple g^ u to matter, and not f^ u , as in the original papers. We note that 
this is not the only possible choice and it would be interesting to explore other options. 



The equations of motion are given by varying the action with respect to g^ u and f^ u : 



R 



2 3 



R 



f^ Y {n)» (yfF*g) +UY(n)» (VT^g 



n=0 
2 ■"• 



-}—T 



(4) 
0, (5) 



n=Q 



where R^ u and R are the Ricci tensor and Ricci scalar due to fav, and 



AC 



The matrices i^^w (-X") are given by 



Y (0) (X) = 1, Y (1) (X)=X-lTrX, 



y (2) (X) = X^ - XTr X + - 1 f (Tr X) - Tr X 



1 



1 



(6) 



(7) 



y (3) (X) = X J -X^Tr X + -X (jTr X) z - Tr X 2 ) - -1 (jTr X) 6 - 3Tr XTr X 2 + 2TV X J J , 

where 1 is the identity matrix. Imposing that T^ v is covariantly conserved, then from eq. Q, the 
Bianchi constraint gives 



0. 



(8) 



n=0 



It can be shown that the corresponding Bianchi constraint given from ^ is equivalent with Q. 
Finally, by performing the constant rescaling 



f, 



M 2 

9 



M 



2 fptv > Pn 
f 



(9) 



we set M 2 to unity. In other words, M 2 was a redundancy that we do not consider a separate free 
parameter. 
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3 Review of background equations 



Cosmological solutions of the Hassan- Rosen bimetric theory were studied in [HI [9], HH [TO! ■ We have 
a separate metric ansatz for each of g^ u and fav (specialized to the case of flat spatial sections): 

ds 2 g = -dt 2 + a 2 (t) dx 2 , (10) 

ds) = -X 2 (t)dt 2 + Y 2 (t)dx 2 . (11) 
The Bianchi constraint given in eq. ^ gives 

^ V + 2^/3 2 + ^/3 3 ) (Y - aX) = 0, (12) 

where overdots denote differentiation with respect to t. There are two options for solving the Bianchi 
identity, which we refer to as Case A and Case B. In Case A, which we will not use, 

Y Y 2 

/3i + 2-/3 2 + ^/? 3 = 0, (not used) (13) 
a gt 

which gives a cosmological solution that is degenerate with GR, for which the fluctuation equations 
reduces to identical equations to those of GR (as noted in [H]). Thus, we focus on Case B, which is 

X = -. (14) 
a 

The Friedmann equations derived from eq. Q and eq. ([5]), together with the Bianchi identity, are 

d\ 2 ( Y Y 2 Y 3 \ 1 

+ m ^ /3o + 3 /3i _ + 3/32 _ +/33 _j =_T °, (15) 

y) 2 + - 2 (ft + ^4 + + = o. (16) 

The acceleration equations can be shown to follow from the two Friedmann equations when using 
the Bianchi constraint. 

In this paper we will only consider the simplest class of solutions, corresponding to /3i = = 0, 
as discussed in [9j . Combining the two Friedmann equations then gives 

tt2 _ P 2 ffofo - 9g 



/3 4 - 3/3 2 3Af| 3 (ft -3ft)' 

y2 g /3q-3/3 2 

a 2 m 2 M| - 3/3 2 ) + /3 4 - 3/3 2 ' 

where H = a/a and p = — Tq corresponds to the pressureless matter density. 



(18) 
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We now define four effective parameters, to be used in our fluctuation analysis, according to 



H. 



ciS 



M 



m 



/3o/3 4 



3/3 2 ) 



3(/3 4 
2m 2 (1 + c 2 ) fo 



Mi 



M: 



A) 



,/34-3/3 2 
-3/3 2 



(19) 



/3 4 - 3/3 2 

The importance of these particular combinations of parameters in the action ^ is as follows. First 
observe in the action that (3q can be thought of as setting the usual g cosmological constant, that /?4 
likewise can be thought of as setting the / cosmological constant, but that the effective "observable" 
cosmological constant that actually appears in (17) is a combination of /3q, /3 2 and ^4. In the 
p — > limit of (17), there is a de Sitter solution, and its Hubble constant H^ s is then related to the 
effective cosmological constant induced by the interaction potential, that is in turned fixed by the f3 n 
parameters. If we also consider p — ¥ in (18), we see that also fn U will have a de Sitter solution with 



possibly different overall normalization, and the parameter c is the proportionality constant between 
g^y and f^ u in the de Sitter spacetime. Further, M 2 is the mass of the spin-2 helicity modes when 
linearizing g^ v and around such proportional background metrics (see appendix [e| . Finally, Mp 
is the effective gravitational coupling constant for p in the cosmological framework (note that this 
will not be the coupling constant for the fluctuations, nor does it necessarily describe the coupling 
in local solutions). 



In terms of these parameters, the H equation can be written 

H 2 



12 9 • 11 i< 



3M p 



(20) 



As usual, the continuity equation for equation of state p = wp with constant equation of state 
parameter w reads 



dlnp 
din a 



-3(1 + w) 



(21) 



with solution p = p$a s ( 1 + w \ Pressureless matter (w = 0) evolves as a 3 , and normalizing the scale 

(22) 



factor at the present time to ao = 1, we can rewrite eqs. (20) and (18) as 

2 



H 



1 



1 



H, 



(IS 



Y 2 



2 (1 + c 2 ) H 2 



+ 



M 2 



H, 



as 



H, 



2{l + c 2 )H 2 s - M 2 ' 



where we also used the relation 



Mp 



M * C 1 + C ) 



M 2 



2H ls 



2(l + c 2 



(23) 



(24) 



terms of the parameters (19). 



Equations (22) and (23) are our final forms of the background equations, expressed entirely in 
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In terms of background cosmology, since the form ( 22 ) is equivalent to ACDM, the only relevant 
parameter is H 2 S /H 2 , which can be constrained by observational data to be close to 0.7, by relating 
it to the usual ratios to critical densities: 

H 2 H 2 

^0 ^0 

It is therefore only the specific combination of f3 n given by the definition of H^s that is constrained 
by the expansion history of the universe, leaving Mp, M 2 and c 2 as unconstrained parameters, 
possibly to be constrained by structure formation data, but there are some further restrictions, as 
we shall see. 

Finally, a comment on the range of these parameters, in particular of the mass parameter M. 
H 2 > Hj s always ( see e.g. fig. pi), we note that if 2 (l + c 2 ) H^ s < M 2 and at some time 



it happens that 2 (l + c 2 ) H 2 (t) > M 2 (which can occur in the early universe), then from (23) the 
/ scale factor Y will be imaginary, which we consider unphysicalj^] Therefore, we demand that 
M 2 < 2 (l + c 2 ) H 2 S . But then, we see that Mp will be negative if also M 2 > 2H 2 S . Negative values 
of Mp would be unphysical, since the matter density would then need to be negative in order to have 
expanding background solutions originating in a hot and dense state (as demanded by observations 
of the cosmic microwave backround). 

To summarize, if we demand that Y should be real and Mp positive, we require 

M 2 < 2Hj s . (26) 

In these bimetric models, we thus need to violate the Higuchi bound M 2 > 2H 2 S [45, 46J already at 
the level of the background. We will comment more on this later, and see also fig. [2j 



4 Perturbations 

We will first consider a general background for g^ u and f^ v . In this paper we will only consider 
scalar gravitational perturbations, except for a brief review of tensor perturbations in appendix [E| 
For the scalar perturbations, following Weinberg [43j we make the ansatz 

ds 2 g = -(1 + E g )dt 2 + 2ad i Fgdx i dt + a 2 ((l+A g )6ij + d i d j Bg)dx i dx j (27) 
ds) = -X 2 {1 + E f )dt 2 + 2XYd i F f dx i dt + Y 2 {{\ + A f )5 ij + didjB f )dx i dx j (28) 

so our set of eight (non-gauge-invariant) scalar gravitational fluctuations is {E g , F g , A g , B g } and 
{Ef, Ff,Af, Bf}. Note that at this point, the ansatz is completely symmetric between the g and / 
metrics, as far as the perturbations are concerned. 

1 It might be interesting to explore this branch of solutions, for example by picking a sufficiently large M that moves 
this region to the very early universe where the current model is in any case not applicable. We will not consider such 
models in this paper. 
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4.1 Gauge invariant variables 



We form gauge invariant combinations of the perturbations in (27). As usual there is no uniqueness 



in the choice of gauge invariant variables (any combination of gauge invariant variables is gauge 
invariant) but we find the following variables convenient: 



2^9 ^ 2 



a 2 B n - 2aF a 



\ E 9 



a 2 B n - 2aF n 



Bt — B„ 



with the definitions 



ci> f 



2 A 



B/ - 2YF 



\ E f 



i 

2A' 



B f - 2YFf 



~x f 



la 2 r 



Y 



(29) 



H = a/a 



K = Y/Y 



(30) 



so is the Hubble function for the / metric. The ansatz ([29j) is roughly speaking "as symmetric as 
possible" between g and /, but complete symmetry is unattainable as the backgrounds are generically 
different. 



4.2 Equations of motion: general background 



Using the gauge invariant variables in the previous section, the equations of motion for the scalar 
perturbations in the g sector become 



1 



m 2 YP 



3 -*/ + ¥ £ 



YK 



F + V 2 B 



1 



2M| 



<5T ° (31) 



(32) 



777 ^ ( 

+ _|p[X($ / -$ fl )-(F^)«] + yQ 



2 + 



+(d 2 + 9 2 )£ 



(33) 



2M„ 2 * 



1 



1 



where j, k are not equal to i. In the /-sector we have 

3 | -¥/ + ¥, 



m 2 a 



^ V 2 ^ + 3^K (ff*, + ¥, j - 



YH 



X 



F + V 2 B 



(34) 



(35) 
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X-^f + K*f)- 2xHaX + Y) 



(36) 



~ £*f + + 3^ + *,) + | 



p 



) */ + ^2 ( 5 J + d D - */) +( 37 ) 

^) + (5| + 5 fc 2 )s]}=o 



with the definitions 



P = (Pia 2 + 2/3 2 aY + /? 3 Y 2 ) , Q = [aft + /3 2 (aX + Y) + f3 3 XY] . 



(38) 



(39) 



At this point we recall that although the separate Einstein-Hilbert actions for g and / in eq. ([!]) 
are of course invariant under separate diffeomorphisms of the two metrics, the mass terms are only 
invariant under diagonal diffeomorphisms that preserve g^ 1 f ■ Thus, the ty g , <E> ff and can 
only appear as differences in the mass terms, and we see this manifestly in the equations. 



4.3 Massless limit 

If we were to turn off the interaction potential in the action, we might expect to find two decoupled 
sets of fluctuations. In fact we observe that if /3i = 02 = @3 = then the two combinations P and 
Q = both vanish, so T and B drop out of the equations entirely, and the g fluctuations and / 
fluctuations constitute decoupled sectors. 

This is a simple observation purely in terms of the fluctuation equations. However, whether the 



fluctuations truly represent decoupled physics is a subtle issue. For example, the condition (14) 
from the Bianchi identity relates the g and / background solutions for arbitrarily small M, but there 
is a priori no reason to impose this condition in the strictly massless theory. (In the language of that 



section, one can revert to Case A, in which case (14) need not be imposed.) But if the condition 



imposed on the background differs between M — > and M = 0, there is a potential "cosmological 
vDVZ discontinuity" [401 141], i.e. the M — > and M = theories could potentially be different no 
matter how small M is taken in the limit. Of course, there could still be a Vainshtein mechanism 
|42] that resolves the discontinuity in the nonlinear regime, but this would not be evident in our 
linear approximation. In figure [5] below, we see some hint of a discontituity, but it is somewhat 
subtle here as we have several parameters to play with. We will not resolve the issue of the existence 
of a smooth limit here, but see also the recent interesting discussions by [28| [29]. 

Now we consider more special backgrounds, first a two-component fluid solution that we will 
refer to as the "exact solution" , then de Sitter and then quasi-de Sitter. 
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5 Exact solution 



It is well-known that in the approximation of a two-component fluid of pressureless matter (dust) 
with equation of state p = (w = 0) and cosmological constant with equation of state p = —p 



(w = —1), the combination of equations ( |22[ ) and (21) admits the exact solution 



3 \ 2/3 

a(t) = ci sinh ( -Hast (40) 



where the constants c\ and -ffdS are 



n ( ^ v /3 = i 1 -^ y /3 > ^ds=y^^ . (4i) 



We can write d = f/a in (15) to express F(i) in terms of a(i): 



= -j^s/H\t)-fcm* ■ a(t) . (42) 



and then using the effective parameters (|19l) we obtain for the scale factor of the / metric: 



2 (l + c 2 )H^ coth 2 H dS t) -M 2 /3 \ 2/3 

Ideally one would now simply use these background scale factors in the fluctuation equations and 
solve them numerically, which would lead to a model for growth of structure in bimetric theory at 
any time t. Unfortunately, we have not been able to complete this program, and instead we will 
focus on special cases and simplifying approximations. 



The simplest special case is that H = constant as in pure de Sitter space, then (43) tells us that 
Y(t) oc a(t), with the constant of proportionality given by 

c=^f. (44) 

We will in general not limit ourselves to pure de Sitter space, but it provides a useful starting point. 
Any departure of the metric from pure de Sitter space breaks the proportionality between the 
g^ u scale factor and the f^ u scale factor. 



6 Pure de Sitter space 

Matter dilutes away as the universe expands, and the universe approaches a de Sitter (pure dark 
energy) solution in the future. To provide some feeling for the numbers, if the evolution would 
proceed according to GR, then it will take around 10 Gyr after present for the exact FLRW scale 
factor a(t) to agree with the de Sitter scale factor aas(t) to within 1% accuracy. 
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For large time, the exact solutions in the previous sections reduce to the approximate solutions 
a(t) -> a dS (t), Y(t) -> y dS (t) where 



«ds(*) = c 2 exp (f/ast) 
*ds(*) = c • c 2 exp (fldS*) 



where 



C-2 



l-Ov\ 



1/3 



(45) 
(46) 



(47) 



and c is the proportionality constant from (43). Although these are of course exact de Sitter solutions 



in their own right, it is useful to consider them as limits of the exact solution for normalization 
purposes. In particular, since there is no Big Bang in pure de Sitter, there would have been no way 
to normalize the scale factor. 



6.1 Gauge invariant variables in dS 



The general gauge invariant variables of (29) have the following dS limits: 

a 2 B 



— A A 



Has 
2 



2aF n 



a 2 B g - 2aF g 



B n 



Defining the linear combinations of field^] 



c 2 $ 



/ 



dS 

2 

a 2 B 



a 2 B f - 2aF f 
f ~ 2aF f 



F f -F a + % 



B 9 ~B 



(48) 



(49) 



f 



(50) 



we are able to separate the scalar gravitational fluctuation equations (31)-(38) into a system of 



massless equations for and massive equations for The T and B fields appear only 

in the massive field equations. The matter perturbations will appear in both sectors. 

6.2 Perturbations in dS 



The equations of motion reduce as follows. In the massless sector: 

1 



r V z *+ + 3H AS H$ + + 



6Tl 
2M* 



(51) 



2 Had we not set M* = 1 by rescaling, it would also enter in these combinations. 
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-$(* + + iT d S* + ) =|§2 (52) 

1 AT 1 * 

*+ + + 3H ds * + + 3H 2 is t + + (d) + 8g) (*+ - *+) = ^ (53) 

We immediately note that these equations are of exactly the same form as the analogous equations 
for perturbations in Einstein gravity. In the massive sector: 

- iv 2 *_ + 3tf ds (ff4>_ + *_) + ^ (i^) (3*_ - ZaH^T + V 2 B) = || (55) 

-A(*_ + ^_) + ?£(I±^ +a *) = § (50) 

+ ff dS (3*_ + 6_) + 3#| S <I>_ + ^ (df + d 2 k ) - *_) + (57) 



^(*--*-)-^(^)^B=^ (58) 



where P and Q are defined in (39). The earlier statement that only differences of the $ and VP fields 
can appear in the mass terms now translates into the statement that mass terms only appear for the 
<&_ and \P_ fields. Thus it must be that the and fields are massless, which is clear above. 
This was also observed in 1251. 



6.3 dS solutions: Massless sector 

We will consider only pressureless matter (dust), both for background and for fluctuations. This 
means that 5Tj = 0. However we will only literally use this condition in this section, since we 
will have sources generating effective pressure and anisotropic stress in the next section. Also, we 
will spatially Fourier transform the perturbations, i.e. write Fourier modes with spatial wave number 
k = |k|. For practical purposes these wave numbers will correspond to wavelengths below the horizon 
scale, so k > H^s- 
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With this understanding, the massless ("plus") gravitational potentials in dS are 

*+,ds = $+,ds = C ie - H ^ + C 2 e- m ^ 1 (59) 

where as before H^s = V^aHq. The matter perturbations are given from the 00 and Oi component 
of the equations of motion: 

= -12^ s C 2 e- 3 ^_ *k 2 (C 1 e-* H ** t + C 2 e- BH <* t ) (60) 
%J = 4C 2 H dS e-^. (61) 

9 

Again these are identical with the corresponding GR solutions. As in GR, there are two integration 
constants C\ and C 2 , and the remaining fields are determined from these. 

6.4 dS solutions: Massive sector 

The massive ("minus") sector can be manipulated to yield a single wave equation from which the 
other fields are determined. Solving for J- and B gives 

2M~ 2 5u + 4H~$_ + 

— ^ ° B < 62 > 

B = ^p^- (63) 

where the anisotropic stress x is defined through 5Tj = didjX- Subtracting the 00 equation from 
the ii equation, and defining 

3 = $_ + _ (64) 
gives the following second-order equation for the massive fluctuation H: 

V 2 3 



with source 



H dS E + ^-f + (M 2 - 2H 2 dS ) E = J (65) 



1 ' - ~" " . l-r2.. TJ 2- 2- 



J= W \ 6P ~ 6(> + 2Hds6u ~ 26il + 3 V X ~ Hdsa * ~ ° *J ' (66) 
where 5p = (1/3)5T£ and 5p = —5Tq. Note that generically, the massive field H is excited by matter 



sources because the definition of 3 in (64) contains g fluctuations that do couple to matter. Because 
of this mixing we cannot think of 3 as purely "new physics" , and in fact it contains a piece that is 
present also in GR, as we will see more explicitly below. 

We note that if we define a new scalar field 

n = ^ (67) 
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its equation of motion is a (sourced) Klein-Gordon equation 

(□ - M 2 ) II = — J . 

where the covariant box operator is 



□ 



d_ 



ou d v 2 

at a 1 



(68) 



(69) 



The identification of the scalar wave equation in massive gravity scenario was previously considered 
in [33], following older work like 



For solving the massive scalar wave equation, it will be convenient to introduce a new "time" 
variable x that goes to zero as t — > oo 



k 



k 



x = = e 

adS-^dS C2-ffdS 



-H dS t 



with C2 from ( 45 ) , and we introduce a rescaled field 

y(x) = yfx H(sc) 



(70) 



(71) 



i.e. y{t) = v /fc/(c2ffds)e- i7dS */ 2 H(t), then y(x) precisely satisfies the inhomogeneous Bessel equation: 

,,2- 



y + -y + i 

x 



2/ = J(x) 



where the ^ parameter (the order of the Bessel function) is 

, 9 M 2 



4 H 2 



(72) 



(73) 



Here M is given in (19) and the source J(x) is found from the massless sector, i.e. from (60) and 



(61 ). Setting x = and Sp = it simplifies to 



jf x) = ^C^sx 3 / 2 | 20c|C 2 ^ s x 3 / 2 ( 2clC2Hl s x 7 ' 2 



k 



k 3 



k 3 



(74) 



We now proceed to write down the solutions of (72) for y(x) and use them to recover II of (68). We 
will only consider M 2 < (9/4)iJ| s here, such that v in (73) is real; the solution for v imaginary is 



discussed somewhat further in the appendix. We give representative plots of the corresponding J v 
Bessel functions in fig. [T] 



6.5 Solutions of massive wave equations in dS: real case 



As we argue in the appendix, the solution of ( 72 ) is 



y(x) = cjJ u (x) + c Y Y u (x) + y p (x) 



(75) 
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Figure 1: The homogeneous massive perturbation S(£) = cj J v {x) j \/x, with x from (701 
and cj a normalization constant, for various v. We note that for v < 1/2 or v imaginary, 
this homogeneous mode diverges for late time. However, for real but small v (the second 
plot) it stays within reasonable values for a several times the age of the universe. We will 
only consider v > 1/2 (as in the first plot) in this paper. 



where the particular solution y p is obtained from (74) and (114) as 



_ 2c 5 CiH dS s 5/2tU (x) 20c 3 5 C 2 H 3 s s 5/2 ^(x) 2c 3 5 C 2 H 3 s s 9/2 ^(x) 
Vp X ' ~ k k 3 k 3 ' 

The s^ )V {x) are special functions that are solutions to the inhomogeneous Bessel equation with 
power-like source, also called Lommel functions, and are defined in appendix [XJ By the asymptotics 
given in that appendix, this y p vanishes as x 7 / 2 for x — > (late times). If we now fix cy = in 



(75), the homogeneous solution J u (x) is proportional to x u for x — > 0. The parameter v depends 



on the mass parameter M and is given in (73). We see that since v < 7/2 (compare fig. [2]), the 
homogeneous solution will be leading (i.e. its leading power of x will be lower) compared to the 
particular solution, at late times. Of course, we can also turn off the homogeneous solution at will 
by setting cj = 0, but the particular solution always remains in the massive scalar IT. 

In fig. [TJ some of the solutions of the homogeneous equation blow up at late time. Moreover, 
this is just the J v solution; the Y v solution blows up at late times for all v in our range. This 
is not necessarily a problem for the theory, though it is a problem for our approximation; all it 
means here is that if these modes would be excited, linear perturbation theory breaks down in the 
(possibly distant) future. Even if nonlinear fluctuations around this background did produce some 
instability in in the future, we do not believe this is relevant to our phenomenological objectives in 
this paper, simply because the far future is more of an auxiliary device here than a region of interest. 
For example, one can picture a scenario where the current model in Hassan-Rosen bimetric theory 
is replaced by another effective theory at late times, where the growing solution is matched to a 
decaying solution in the new theory. It would be interesting to learn that this is not possible and 
actually the presence of these modes rule out the theory for some values of M, but at the moment 
this is not clear to us. See also the conclusions for more comments on this. 
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We summarize in figure [2]why v > 1/2 seems the only reasonable choice in the current model. 



M 2 P < 



one — > oo both — > oo 



M 2 /-ffjs ■ 5/4 2 9/4 2(1 + c 2 ) 



^ = 3/2 z, = 1 i/ = 1/2 zy = 



Higuchi v i ma g inar y Y imaginary 

bound 

Figure 2: Summary of some interesting values of the mass parameter M and the corre- 



sponding v values from ( 73 ) . We indicate when only one homogeneous mode goes to infii4ity 
at late time, or when both modes do, and when v becomes imaginary (which is not a prob- 
lem in itself, but see fig [TJl and when Y and Mj> become negative (which are big problems 
for the background in the current model). 



6.6 Parameter values 

We now have four integration constants that need to be specified: First, the constants G\ and C2 



of (59) for \l/_|_ i( jS) which we can think of as initial conditions (ICs) for ^t+^s an d VP+dS) i- e - the 
gravitational potential perturbation of the massless sector. Or equivalently, we can think of them 
as ICs on the values (but not the derivatives) of (Jpds and Su^s- Second, we have the the constants 
cj and cy from the homogeneous solution for Eds, which we can think of as ICs on Sds and Hds, i-e. 
the "wave-like" field in the massive sector. Thus in general we have four parameters {Ci, C2, cj, cy}- 

For convenience we give the mass M 2 and the comoving wavenumber k 2 in terms of the Hubble 
parameter value H^s = V^aHq that is asymptotically approached in the future, i.e. not in terms 
of the Hubble parameter value Hq of today. 



6.7 Conclusions in dS: nothing new for cosmological geometrical probes 



Because the perturbation equations for the "+" subscript fields, eqs. (51)-(54), are identical to the 
GR equations, and the "— " subscript fields only appear in the massive equations (55 )-(58 ), and never 
appear in the GR-like equations at this order, the matter perturbations and hence the linear growth 
of structure will be identical to that seen in general relativity. It is important for this conclusion 
that the matter perturbations are completely determined by the 00 and Oi components of Einstein's 
equation, which are of course constraints. In general in massive gravity this may not be the case; 
see the conclusions for related comments. 

We also note that the argument in the previous paragraph leaves much to be desired from a 
phenomenological standpoint. For example, as light rays propagate along geodesies of the g metric, 
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presumably one cannot form linear combinations of fluctuations in the electromagnetic field as 
we have done here that look identical to the GR equations, even in pure de Sitter, so one could 
expect that bending of light and similar local experiments will be grossly affected even where the 
cosmological matter perturbations are not. However, we emphasize that local solutions are not yet 
well understood, and we follow the practice of [9] of lumping these unknowns into Mp, that is never 
directly relevant in global equations. 



7 Quasi-de-Sitter space 



We have seen that bimetric cosmological gravitational perturbations in de Sitter space are identical 
to those of GR. For this and other reasons, it is of interest to consider a universe that represents 
a small devation from from the de Sitter background, and then consider cosmological perturbation 
theory in this slightly generalized background. This kind of approach works well in inflation, but it 
is not often used for late-time cosmology where the presence of matter complicates things. We will 
see some of these complications, and why the quasi-de-Sitter approach is still useful for our purposes. 
For some background on quasi-de-Sitter, it is useful to consult a review of inflation, e.g. |6|. 



7.1 Background 

We define quasi-de-Sitter space as near-exponential expansion of the scale factor 

e = -~>0 (77) 

where as usual H(t) is defined as H = a/a, and the square in the denominator makes e dimensionless. 
For pure de Sitter space, H is strictly constant so e = in pure dS. (The reason for the tilde will 
become apparent shortly.) We remind the reader that also in inflation e is often first defined in 



terms of the geometry, just like in (77). The inflationary slow-roll parameters in inflation are then 



given in terms of some scalar field potential, for which we have no direct analog here. 
It is a simple matter to compute —H/H 2 from the Friedmann equation (22) to obtain 

3 

6 ~ 2Q A a 3 



(78) 



where we have used the zeroth-order relation H(t) ~ -f^ds- We see that e encodes the fraction of 
matter in a universe dominated by dark energy. Clearly, this should be a small parameter in the 
future, and as we approach matter domination, the quasi de Sitter approximation breaks down in 
the past. 



3 In fact, the expansion in e.g. Ch. 8 of Weinberg [43] using the correction factor C(x) is the inverse expansion of 
this, there it is the ratio of dark energy to matter, which breaks down around present and more severely in the future. 
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We split the Robertson-Walker scale factors of the two metrics into products of de Sitter scale 
factors and correction factors: 



a(t) = q a (t)a dS (t) , 9o (i-K») = l (79) 
Y(t) = q Y (t)Y dS (t), q Y (t -+ oo) = 1 (80) 

so that the functions q a (t) and qy(t) capture the "quasi- ness" of the expansion. Let us consider the 
exact solution (40), (45) for the two-component fluid, then we can extract the quasi- ness for a(t) as 

\ 2 / 3 / i \ 2/3 



q a {t) 



a{t) 



«ds(*) 

exactly, where now we define 



2 2/3 sinh (l g dst)' 
exp H dS t 



1 



e(t) 



e(t) = 6e 



-3H dS t 



(81) 



(82) 



which agrees with (78) to lowest order (hence the tilde in (|78|)). For some numbers to keep in mind, 
e ~ 0.01 at Hot ~ 2.5, e ~ 0.1 at H^t ~ 1.6 and e ~ 0.5 at H^t ~ 1. We will prefer to stay at 
e < 1, which limits us to Hot > 0.7 as a matter of principle. (In actual examples, we will find greater 
limitations than this.) For small e, we can expand the quasi-ness of a(t) in e: 



q a (t) 



l 



1 , x 1 



9 w 324 

We can then easily compute the (square of the) Hubble function: 



a 2 



H, 



(IS 



2 2 9 
1 + 3 e + i £ 



leading to 



H 



<*) + 



(83) 



(84) 



(85) 



Note that because we defined e(t) as (82), this is not exactly e(t) of (77). This distinction is one 
of convenience and merely amounts to a rearrangement of higher-order perturbation theory in the 
"true" quasi-de Sitter parameter e. 



For the / metric one can argue similarly. From (43) and with the dS solution in (45), 



Y(t) 



Ids 1 



9 



^6 2 + ... 
324 



where 



Cy 



cy,2 



M 2 + 4(1 + c 2 )Hj s 
M 2 - 2(1 + c 2 )H 2 



(86) 



(87) 



(IS 



M 4 + 44(1 + c 2 )M 2 H 2 s - 20(1 + c 2 ) 2 Hj s 



(M^ - 2(1 + c 2 )H%) 2 

This qdS expansion of the scale factor Y captures the loss of proportionality between Y(t) and a(t) 
as we leave the pure dS regime and enter the quasi-de Sitter regime. (Note that the constant cy 
is never unity.) We summarize the results for qdS expansion coefficients for the various derived 
background quantities in appendix Oat linear order, which is all we will use explicitly in this paper. 
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Figure 3: Background solutions for the Hubble functions H = a/a and K = Y/Y, with the 
latter for c = 6 and c = 1/6, cf. (44). (Here the c = 1/6 curves are included for illustration 



only, we never use values for c this low). Linear (in e) qdS in black, quadratic qdS in dashed 
red, exact solution in dotted blue, and H t ~ 1 is roughly present. We see that for times 
Hot > 0.7, qdS remains a good approximation for the g background, and also for the / 
background for c = 6. 



7.2 Perturbations in qdS 



We write the general perturbation equations in section |4.2| as D mn <p n = J m for a collection of fields 
(f) m enumerated by m = 1 , . . . , 8 and a differential operator D mn and sources J m . We organize the 
expansion as follows: 



D — /5 U 4-fD e 



(89) 



from which we write the first order equation as 



e(i7m 



+ 7 e 



»0 



(90) 
(91) 



using the zeroth order de Sitter equations D^n^o = Jm, and neglecting quadratic order in e. As 
expected, the zeroth order (pure de Sitter) fields act as additional sources J m for the first order 
quasi-de Sitter fields. It is useful that from this vantage point, the differential operator on the 



left-hand side of (91) is the unperturbed de Sitter differential operator Df : 



Note that (91) contains terms with the time derivative e of our perturbation parameter. In 



principle by using (82) we can express these entirely in terms of e, and the latter would then drop 



out of the equation. In practice, it is convenient to work with the unperturbed dS differential 
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operator also in qdS, so we will not substitute in for e and instead simply write 

<T = fife + <ds • (92) 

The price to pay for this convenience is that viewed as expansions in x ~ e~ Hdst , it is not guaranteed 
that every term in <?^g is suppressed compared to every term in </>^g, and in fact it will generically 
not be the case, but the overall series in x does display relative suppression. We summarize this fact 
m table [l] For more explicit comments on this issue in the simpler setting of ordinary GR, we refer 
the reader to appendix |D] 

It would be useful to know when the linear (in e) qdS approximation is valid to some some 
prescribed accuracy for the perturbations. However, we currently do not have perturbations in 
an exact solution to compare to in bimetric theory, so it is hard to be absolutely precise about 
this (and if we did have an exact solution, the question would be rather pointless). As a first 
check, we compute the relative errors of the linear and quadratic qdS approximations versus the 
exact solution in cosmological perturbation theory in pure Einstein gravity in the aforementioned 
appendix [D] As a second check, we have performed some preliminary analyses of the quadratic (in 
e) qdS approximation also in bimetric theory, in particular how the bimetric perturbations in the 
quadratic qdS approximation differ from the linear qdS approximation (typically if second order 
perturbation theory produces significant changes, perturbation theory has broken down). To be 
clear, for the purposes of this paper we only use the quadratic qdS approximation for auxiliary 
checks and we do not display it in plots. We find numerically that the linear approximation in e is 
good to about 10 to 30% for the bimetric perturbations (depending on the field) for H^t > 1.3 — 1.5. 
This will be indicated by shading the region below this in the plots. 

From the good accuracy of the qdS background in section [3] one could have hoped that the 
linear qdS perturbations would have extended further back than H^t > 1.3 — 1.5, since the future 
is of no direct use for phenomenology, but there was of course no guarantee that this would be 
the case. On the good side, since going to quadratic order seems to give some improvement in 
our preliminary analyses, we believe that the qdS approximation at higher orders will be useful 
also for phenomenological purposes, and not only as indirect checks of numerical solutions of the 
perturbation equations in the exact background. 



7.3 Parameters 

We use wavenumbers k = WH^s, k = (5/2)Hds arid k = (l/2)#ds as representative cases, the latter 
only for internal checking of the analytics, as we will describe later. For k = (5/2)H^s we will impose 
the following values on the matter perturbations at H$t' = 1.5. (One would have liked to do this 
at present t = to, but the qdS approximation needs to be valid in the region where we set initial 
conditions.) We obtain the values from GR (see appendix |P]) : 

^ = 8.93 • 10- 5 , ^ = 0.72 • 10~ 5 (93) 
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For k = (l/2)HdS> we have instead 

^ = 2.43 • 10" 5 , ^ = 0.72 • 10- 5 . (94) 
M 2 M 2 V ' 

We view these as having roughly 10% accuracy. It is nontrivial to extract the values directly from 
data, but doing so would be useful in a more phenomenological analysis, rather than comparing 
directly to Einstein gravity, since we are of course modifying gravity. 



7.4 Analytical series solution 

The qdS equations are more complicated than the dS equations, but as the differential operator on 



the left-hand side of (91) is the unperturbed de Sitter differential operator, the general strategy for 
solving the differential equations is the same. In particular we again find a massive inhomogenous 
Bessel equation for H q ds, just with more complicated sources. We will not turn on this homogeneous 
solution, since if we did, it should have been included at dS order (if this is unclear, it may help 
to consult our analogous comments in GR in appendix |d| . We introduce y = x 1 / 2 , cf. (70). We 



compute series expansions of the right hand sides of all the perturbation equations to see which 
powers of y actually occur, and arrive at the following series ansatz: 

= $( 7 )y 7 + $(8)y 8 + $(ii)y n + $(i2)y 12 (95) 
= ^(7)y 7 + ^(8)y 8 + ^(ii)y n + ^(i2)y 12 (96) 

dp = 5p(j) y 7 + 5p (8 ) y 8 + 5p ( ii) y 11 + Sp {12) y 12 (97) 

Su = 5u(j) y 7 + <kt( 8 ) y 8 + Su^ y 11 + fai (12 ) y 12 (98) 

5 = E (7) y 7 + E (8) y 8 + E (n) y u + E (12) y 12 (99) 

$_ = $_ (7) y 7 + $_ (8) y 8 + $_ ( ii ) y 11 + $_ ( i 2 ) y 12 • (100) 

(The J- and B fields are determined from these, as before.) The explicit expressions for the coefficients 
obtained in this way are not terribly illuminating so we do not present them in full, but to give an 
idea of what they look like for our fixed parameter values (in particular v = 1), we find coefficients 
of the rather manageable form 

_ 6c 2 ff d 3 s 48M 4 - 176M 2 ff d 2 s + 205ff d 4 s 128^ 
^ • (M 2 -74F d 2 s )(M 2 s -2F 2 s ) " CJ 97 P ' ° J ^ 

* 8c|ff d 4 s M 2 + 34ff 2 s = _376c|^ 

(8) 3fc 4 (M 2 -74tf 2 s )(M 2 s -2tf 2 s ) 1 291 k 4 1 { ' 

and so on. For our parameters, 2(1 + c 2 ) = 74, so both denominators are (M 2 — 2(1 + c 2 )H 2 s )(M 2 — 
2ffds)i i- e - the expansion breaks down not only for early times but also for two of the distinguished 
mass parameter values in figure [2I The breakdown point M 2 = 2(1 + c 2 )H^ s could have been 



anticipated from the background expansion (87). 



For k = (l/2)Has, we find that the present time t = to corresponds to x = 0.47 (see (70)) 



and for k = (5/2)Hds we find x = 2.37. So for the larger k values, the more phenomenologically 
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Table 1: Leading expansion powers in the analytical solution. The notation means that 
the leading term at late time is x A . 



interesting ones, we expect this analytic version of the qdS expansion to break down, even when a 
numerical qdS approach would still be valid. (However, the analytics may be somewhat better than 
this, since many terms are actually expansions in x/k, which is independent of k.) Therefore we 
focus on k = (1/2) H^s when we use the analytical series solution. In all, the lessons we learn from 
the analytical solution is that there are certain degenerate special parameters, and we can quantify 
when a given approximation breaks down at least for small wavenumber k. None of this will be 
evident in the following, since we have already identified useful parameter values and we will not 
bother to show plots comparing the analytical and numerical results, we will just state here that 
they agree to the extent we expect them to. Perhaps the most important use is as cross-check with 
the numerics for low k. We now turn to the numerics. 



7.5 Numerical solution: general 

Naively we would expect that the energy density would receive a slight positive correction since 
going away from pure de Sitter expansion means that there is less expansion and the friction term 
due to the expansion is thus smaller. But in bimetric gravity the situation is more involved since 
both the massless and massive background sector will contribute to the correction. For example, the 
initial conditions set for the massive wave (i.e. the two integration constants for the homogeneous 
solution) will affect the correction to the energy density. 

One way to fix parameters would be to use p, u and p at present to fix ICs Ci, Ci and cj (we 
always set cy = 0). Another way is to fix cj = and fix C\, C2 from p, u at present, which is what 
we will do here (with the exception of figure [7]) . For the values (94) for Sp and 5u, and k = (5/2)Has, 
we obtain 

Ci = -2.19 • 1(T 5 , C 2 = 3.33 • 1(T 5 , cj = . (103) 
where we have normalized the gravitational potential as in GR (see appendix [Pj). 

7.6 Numerical solution: gravitational potential 

We first show two plots of the gravitational potential in figure [4} One generally expect that for larger 
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wavenumber k, the natural time variable x in (70) is larger, so if we were to series expand the Bessel 
(and Lommel) functions, we would need to keep more terms. In other words, for large wavenumber 
k, the solutions "explore" the Bessel functions more, and the oscillations there can carry over to the 
gravitational potential. 



7.7 Numerical solution: density contrast 

We form the density contrast 

8 =. 6 P-* H6u (104) 
P 

where p is the background matter density. This is used to compute the growth factor. 

We observe that although the field depends on M, the density contrast 5 seems to depend 
on M much less. This is partially because of the way we fix boundary conditions, which is imposed 
directly on 5 and therefore only indirectly on \]/ + . Nevertheless, although the 5 we see here does 
not differ appreciably from GR around H^t > 1.3, it does differ in the future, so we would expect 
that if we go away some time from the point at which we give the ICs (here Hot = 1.5), there would 
in fact be some controllable difference, which is where the phenomenology of matter perturbations 
could begin. 

It is of great importance whether there is a Vainshtein-like mecanism here. If there is a finite 
gap between the GR solution and the bimetric solution for any value of M, one might be tempted to 
conclude that there is in fact such a mechanism at work. However, there are also other parameters, 
for example cj, that can be turned on to try to mimic GR at zeroth order. See figure [7| 
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Figure 5: Numerical qdS plots of the density contrast 6 for various wavenumbers k. Left 
panel: small k values, k/H^s — {5/2, 3, 7/2, 4}. Right panel: intermediate k values, k/H^s = 
{19/2,10,21/2,11}. Each curve is plotted for two M values, M 2 /Hj s = {1/50,3/2}, but 
the curves for the two M values are nearly coincident in some cases. The shaded area is our 
estimate for when the qdS approximation breaks down. 



7.8 Numerical solution: massive wave 

The massive wave LT/S is plotted in Fig. [6j Again, the "bumps" are in the region where the 
approximation has already broken down, and so cannot be trusted. Still, also here we learn something 
about the massive wave around H$t ~ 1.5, and we see that it does depend on M, as one would expect. 

One can use the existence of the homogenous mode to see if one can recreate GR. We show some 
simple attempts in this direction in figure [7} 

8 Towards ACDM 

We see that it is in principle possible to obtain good accuracy with the quasi-de Sitter approximation, 
but what suffices depends on the detailed application. In particular, with our current understand- 
ing it seems that it would be beneficial to automate an arbitrary-order qdS approach in symbolic 
manipulation software, especially if one wants to go to relevent eras such as z ~ 1. This is certainly 
possible, but outside the scope of this work. Even better would be if one could solve the fluctuation 
equations in the exact solution. At the moment we do not know if this is feasible in practice. 
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Figure 6: Numerical qdS plots for the massive scalar gravitational fluctuation IT = ^/a 2 , 
for wavenumber k — (5/2)iJds (left panel) and k = 10-f/ds (right panel), and for mass 
parameter M 2 /H 2 S = {1/50,4/5,5/4,3/2}. 
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Figure 7: Numerical qdS plots for the density contrast S, for the homogeneous mode in H 
turned on, i.e. cj 0, unlike in the other plots in this section, where cj = 0. 
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9 Conclusions and outlook 



We have computed the scalar fluctuation equations in de Sitter space and quasi de Sitter space, and 
we have found some analytical and some numerical solutions. There is much left to do as regards 
phenomenology. 

It would also be very interesting to perform the analogous Hamiltonian analysis to analyze linear 
and nonlinear stability of these fluctuations. In GR, we can solve the gravitational perturbation 
separately, which then completely determine the matter perturbations. It is not clear that this 
strategy should automatically work here, as some of the constraints may become dynamical. But 
at least in quasi-de Sitter, there seems to be no real issue with this. To completely understand this 
issue, we would also need to perform a Hamiltonian analysis, which is beyond the scope of this work. 
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A The inhomogenous Bessel equation 
A.l General 

Let us begin by recalling some elementary facts about general inhomogenous 2nd order ODEs to set 
notation: 



for some polynomials p(x), q(x). Let y\ and 1/2 be fundamental (linearly independent and normalized) 
solutions to the homogenous ODE D[y] = 0. We set for the general solution 



for two unknown coefficient functions u±, U2- Using variation of parameters we find solutions for 
u\(x) and U2(x): 



D[y] = y" +p(x)y' + q(x)y = g{x) . 



(105) 



y g (x) = ui(x)yi(x) + u 2 {x)y 2 {x) 



(106) 




(107) 



where the Wronskian W is the usual determinant 



W[yi,y 2 ] = yiy 2 - yiy'\ ■ 



(108) 
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A. 2 Bessel 



For Bessel functions, the fundamental solutions are yi(x) = J u {x) and 1/2 0*0 = Y u (x) and the 
Wronskian is quite simple: 

2 



W[J V ,Y V ] = J V+1 Y V -J V Y V+1 



7TX 



so the coefficient functions become 



ui(x) = —— dxY u (x)g(x) ■ x , u 2 (x) = - / dx J u (x)g(x) ■ x 



(109) 



(110) 



where again g(x) is the right-hand side of the inhomogenous equation, and the general (and generic) 



solution is simply (|106|): 



ui{x)yi{x) + U2{x)y2{x) 

~ I dxY u (x)g(x) ■ x J J u (x) + ( ~ / J v (x)g(x) ■ x ) Y„(x) . 



(Ill) 
(112) 



A simple way to specify ICs is to fix in (112) and write 

V = Vh + V s , (113) 

with the usual two free integration constants in the homogeneous piece y^, and no free constants in 
y g . Then (112) above represents the particular solution. If we specialize to a power-like right-hand 



side, g(x) = x M , and fix xq = 0, then 

^/particular \x) — 

where s„ v {x) is a Lommel function. The series expansion of this Lommel function at x = is 



| £dxY u (x)x» +1 ^j J u {x) + (| £dxJ v (x)x" +1 ^ Y v {x) (114) 

(115) 



^/i,v(-^) 



X' 



+ 0(x 



/H-3\ 



(116) 



( At _ I/+ l)(^ + l/+ l) 

i.e. the order of vanishing at x = is independent of u, which is not evident from the integral 



representation (114) 



B Massless limit 



In the limit M -> (0 2 -> 0) we find 



^3/2(2;) 



2 sin x — x cos x 

372 » y 3/ 2 (^) 



x" 



2 x sin x + cos x 

7T X 



3/2 



(117) 



Using this, we find from (114) that the associated Lommel functions reduce to simple powers: 

x 2 + 2 

S3/2,3/2(a?) = Zzlo ( 118 ) 



•55/2,3/2 (" 



X 3 / 2 
,3/2 



(119) 
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C Replacement rules 



The range of t is such that the dimensionless e(t) is considered small. In the linear e expansion, we 
have for the scale factors 

a(t) 
Y(t) 

with 

a € = 
c Y = 

and for the derived quantities 

H(t) EE 
K(t) EE 
X(t) EE 

where the coefficients are constants 

H e = 3a e Y e (127) 
if e = c Y /3 (128) 
X e = 2a e {Y e - 1) . (129) 

These simple expressions are sufficient to eliminate all derivatives on the background. 

D Quasi-de-Sitter expansion in Einstein gravity 

In this appendix we apply the qdS expansion to scalar perturbations in Einstein gravity. We compare 
the results of first and second order qdS expansions to the exact solution for the two-component 
fluid (with only dust and dark energy, which should be a good approximation to physical cosmology 
in this time interval). 

The well known GR scalar perturbation equations for dust with the conventions used in this 
paper can be written as: 



= a dS (l-a e e) 
= *ds(l - a e c Y e) 



(120) 
(121) 



1 
9 

M 2 + 4(1 + c 2 )H% 
M 2 - 2(1 + c 2 )Hj s ' 



(122) 
(123) 



e - = H dS (l + H e e(t)) 
e ^=H dS (l + K e e(t)) 



Y 



c(l + X e e{t)) 



(124) 
(125) 
(126) 
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* + H§ + 3H ( H$ + if) + 2H$ = 
1 



2a 2 



d l dj ($-*) = 



Recall from eq. (83) that the scale factor expanded to second order can be written as: 

a(t) = a dS (t) (l - ~e(t) - ^e(tf 

Using this expansion, the equations of motion for the perturbations can easily be recast in the general 
form 

Dirf = (p% + Dfj + Z>g) (4 + 4 + <jP £2 ) =J° + JF + jf 
which we split order by order into three sets of equations: 

D%4 = J? 

d%4 = Ji- Dfrfl 

It is possible to analytically solve these equations, obtaining for the gravitational potential: 

^ = Cie- Hdst + C 2 e- 3Hdst 

y £ = A (5C ie 2H ^ + 3C 2 ) e~ 6 ^* 
15 

$ £2 = _ (50Cie 2HdS * + 27C 2 ) e - gHdst 
90 

where we have set the integration constants of the 1st and 2nd order equations to zero. In fact, 
the solutions of the homogenous versions of these equations are effectively of zeroth order, so to 
be consistent with the qdS expansion we should turn them off. With this understanding, the three 
expressions above are nicely separated in order. In the language of the main text, we have A<j, = 1, 
A«j, £ = 4, A^ 2 = 7 (which is a compact way of stating that the leading terms for large t are e~ Hdst , 
e -4H AS t anc j e -7H dS t^ reS p ec tively) . Here the suppressed terms in each fluctuation are less suppressed 
than the next order, as one would expect. This will be different for the matter perturbations below. 



It is then straightforward to compute energy density and velocity perturbations using (130) and 



(131). For the energy density we obtain: 



M 2 15 



4 ^35^| s C 1 e 2 ^_^ C 7 1 + 225F2 s C 2 -^C 2 e- 2 ^A e -6H dst 

c 2 C 2 / 
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5p £ 



M 2 p 



45 



7560tfJ s Cie 2HdS * 



C 2 e 



We observe that A<5 po = 3, A$ Pe = 4, A$ p 2 = 7. Because the fields begin to mix at first order 
in the qdS approximations, also fields that are suppressed at zeroth order, as Sp is, receive a first 
correction that is relatively big if the other fields are relatively big. In particular, 5po starts at 
e -3H AS t an( j linear qdS field has a e~ 4Hdst piece. Moreover, there can be terms in the lower order 
fields, just from solving the equations, that strictly belong to higher orders in the expansion; this is 
the case for the e~ 5Hdst term in 5p$. We will typically keep such terms at the order at which they 
appear, but they cannot be considered reliable for truncation at the given order. 



For the velocity perturbation: 



Sup 
M 2 p 



5u £ 

mJ 



5u f 



Ml 



AC 2 H dS e- 3Hdst 
12H dS (Cie" 4 ^ 3 * + C 2 e- 6Hdst ) 
-#ds (40Cie~ 7HdS * + 29C 2 e~ 9Hdst ) 



for which Ag uo = 3, Ag Ue = 4 and A$ u 2 = 7, like for bp. 

To compare these results with the ACDM solutions, in the following called fr/, 5pf, and Suf, we 
have chosen C\ and C 2 such that fr f asymptotically matches fro in the future. Notice that this is a 
slightly different choice of integration constants from that used for bimetric theory in the main text, 
but it is more suited to this analysis. Typically the two choices give very similar results. 

With the initial conditions fr/(i*) = 1CT 5 , */(**) = , where i* is HqU = 0.002 (during 
recombination) we found 4 ' 



C\ ~ -2.28 • 10" 



C 2 ~ 4.89 • 10" 



(132) 



Perhaps the best way to get an intuitive idea about how good our first and second order qdS 
approximations are is to consider ^ plots like those of figure |8| A more precise measure is the 
relative error: 



/ 



~ 0.015 , 



H t=l 



/ 



(* + + 



0.009 . 



H t=l 



so both the first and second order expansions are good to about 1% around present. For our purposes, 
it is also important to have an idea when the approximations break down. We find that 



(fro + % 



0.1 



fr/ - (fr + ^€ + * e 2 

fr/ 



0.1 



H t=0.5 

0.5 for the first and second order qdS 



H t=0.7 

i.e. we are down to 10% accuracy at Hot ~ 0.7 and H^t 
approximations, respectively. 

4 We have not been careful with the overall factor here, since we do no actual phenomenology in this paper. If 
the factor changes, all fields would simple be multiplied by the same correction factor, since we are doing linear 
perturbation theory. 
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Figure 8: Comparing approximations for the GR gravitational potential ^. The first order 
qdS approximation is black solid, the second order qdS approximation is red dashed, and 
the exact ACDM solution is blue dotted. 



The range of validity of the various approximations for the matter perturbations might vary with 
k. We analyzed what happens if 1/k 2 lies between the horizon scale and two orders of magnitude 
below the horizon scale, i.e. when: 

H 2 S < k 2 < W0H 2 s 

Around present, for all k in this range, the approximations do not differ more than 10% from the 
exact solution, i.e. 



0.1 < 



Spf - (dpo + Sp e ) 



Sp f 



< 0.01 , 0.007 < 



H t=l 



Spf - (Spo + 5p e + Sp e 



Sp f 



< 0.005 . 



H t=l 



To be more precise, we found a small range of k in wich the density perturbation qdS expansions 
are as good as the \E f expansions or even better, see figure [9| and |10[ in particular we find: 



3#ds < k 2 < !2H dS 
H t > 0.7 



Sp f -(5po+Sp t ) 



< 0.1 



3#ds < k 2 < 12H dS 
H Q t > 0.5 



8pf-[Sp +8p c +Sp e2 ] 



< 0.1 



Having computed density and velocity perturbations qdS expansions, it is simple to compute the 
corresponding expansions for the comoving density contrast 5 = (5p — 3H5u)/p: 

5i = <5o + 5e 62 = 5q + 5 e + 5 e 2 
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Spj-(Sp +Sp e 
S Pf 0.3 n 



6pf — {5pa+t>Pe+6p c 2 ) 
6p f „, 





0.6-^ O.J '•"6.8 0.9 1.0 



Figure 9: Signed percentage differences between density perturbation expansions and exact 
solutions for various k values. Left panel: differences between 1st order qdS and exact 
solutions. Right panel: differences between 2nd order qdS and exact solutions. In the 
interval 3i?as ^ k 2 < 12i?ds (solid blue lines) the percentage differences are always smaller 
than 10% back to Hot ~ 0.7 in the case of 1st order qdS, and back to Hot ~ 0.5 in the case 
of 2nd order qdS. 



Concerning the validity of 6 approximations we find the same results as for 5p above. This is simply 
because the errors in the 6u expansion are smaller than the errors in the bp expansions for all values 
of k 2 . 

Concerning the growth index we find that in general, the values computed with our 1st and 2nd 
order qdS approximations are trustable only for Hot > 1.5 It is not hard to understad why this 
happens. We recall the explicit expression for the growth index: 



7 



In/ 



In (O m /a 3 



with 



din 5 
din a 



Now, given errors in 5 and 5, it is straightforward to compute the expected error in 7: 



A7 




ln(tt m /a 3 ) 



In general for H^t < 1.5 we have |A<5| ^> \A5\. This is why the growth index expansions start 
deviating from the exact function before (in the sense of coming from the future) the other quantities 
of interest, i.e. the approximation for 7 itself is a little worse than that for e.g. ^. 
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Figure 10: Density perturbations for k 2 = 5-ffds and velocity perturbations: the first order 
qdS approximation is black solid, the second order qdS approximation is red dashed, and 
the exact ACDM solution is blue dotted 

E Bimetric tensor modes in de Sitter 

We are not directly interested in tensor models in this paper, but we would like to compare our mass 
parameters to the mass parameter of the tensor fluctuation. Consider gravitational waves traveling 
in the z-direction. The conditions for tracelessness and divergence-freeness of the perturbations are 
solved by the following ansatz 

ds 2 = -dt 2 + a{tf {dx 2 + 2h 9 xy (t, z)dxdy + h 9 xx {t, z)(dx 2 - dy 2 )) , (133) 
ds 2 = -c 2 dt 2 + c 2 a(t) 2 (dx 2 + 2h f xy (t, z)dxdy + h f xx (t, z)(dx 2 - dy 2 )) , (134) 

where a(t) = C2e Hdst . In terms of the linear combinations 

^xx = hxx ~l~ c h X x > hxx = ^xx ~ hxx > (1^5) 

^Xy C , h x y ~ ^y , (l36) 

the linearized equations of motion become 

d 2 d 1 d 2 \ 

W + ZH dS - - - 2 — 2 + M 2 j h XXtXy = , (137) 
d 2 d 1 d 2 \ 

w +3HdS dt~a^d^ h+ ^y = ' (138) 



where 



1 



M 2 = [l + ^)m 2 (eft + 2c 2 ft + c 3 /3 3 ) , (139) 
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Figure 11: Comoving density contrast 5 and growth index 7 for k 2 — 5ifds : the second 
order qdS approximation, dashed, the first order, solid, and the A CDM solution, dotted. 



which can be compared to (19), and where we used the background g and / equations to reexpress 
00 and j3i as 



c 4 /3. 



4 



3H 2 /m 2 - 3/3ic - 3/3 2 c 2 
3c 2 H 2 /m 2 - fac - 3/3 2 c 2 



/3 3 c 3 
- 3/3 3 c 3 



(140) 
(141) 
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